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Abstract 


System of linear equations plays an important role in science and engineering. One of the applications of this system 
occurs in the discretization of the partial differential equations. This paper aims to investigate an experimental 
compatison between two kinds of iterative models for solving the elliptic partial differential equations. Different tools 
of solution such as stationary and non-stationary iterative methods with preconditioning models have been studied. Two 
types of discretization schemes (centered and hybrid) have been also considered for the comparison of the solution. 
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Introduction 


Solving of system of linear equations (SLE) often occur in a wide variety of area, including numerical 
differential equations, eigenvalue problems, economics models, design and computer analysis of 
circuits, power system networks, chemical engineering processes, physical and biological sciences [1]- 
[4]. Such systems are typically solved by two different types of methods; iterative methods and direct 
methods. The nature of the problem at hand determines which method is more suitable. 


A direct method for solving large sparse linear systems of the form of Ax=b involves explicit 
factorization of the coefficient matrix A into the product of two certain matrices such as triangular 
matrices [5], [6]. This is a highly time and memory consuming step; nevertheless, direct methods are 
important because of their generality and robustness. For linear systems arising in certain applications, 
they are the only feasible solution methods. 


Corresponding Author: s.a. edalatpanah@gmail.com 


4^5 https: //doi.org /10.22105 /cand.2022.155122 


Furthermore, direct methods provide an effective means for solving multiple systems with the same 
coefficient matrix and different right-hand side vectors because the factorizations need to be performed 
only ones. The trouble with sparse direct methods is that, because fill-in, their memory requirement 
grows faster than the size of the problems [7]. 


Iterative methods are generally much simpler to describe in details than direct methods, because the 
lack of fill-in and simple matrix access remove the need for sophisticated data structures and graph 
theoretic analysis. On the other hand, the numerical behavior of iterative methods is much more 
complicated than that of direct methods [7]-[10]. 

To sum up, in sparse iterative methods, the spectral properties of the coefficient matrix determine the 
effectiveness of the method. Therefore, it is a challenge to make them robust enough to serve in portable 
libraries and environment used for a wide variety of problem domains. In this paper we apply two 


preconditioned iterative method for solving PDE problems and also compare these two kinds of 
strategys. 


2 | Iterative Methods for Systems of Linear Equations 


Consider system of linear equations as 
Ax =b, (1) 


A system of nonlinear equations, given in fixed point form 
x = T(x), (2) 


can be solved by choosing an initial approximation x € IR" and by successively calculating more 
accurate approximations using the iteration: 


x eT. Rud s (3) 
Stationary iterative methods for the Eg. (1) have the basic form 


x(9 - Bx€-D 4c, k 2 125,.., 3) 


Where B € R"*" and c € IR"do not depend on k . Initially, however, the linear system Eg. (1) must be 
transformed into a form to which the iterative method is applicable. For instance, the coefficient matrix 


can be written as: 
A=D-L-U, (4) 

where D = diag(ay, ... , Ann), L is a strictly lower triangular matrix and U is the respective upper triangular 

part of A. For any splitting, A=M-N, where M is nonsingular, the iterative method for solving linear 


systems of Eq. (7) is 


x) = MINED + M-1b, k 2125, .., (5) 


This iterative process converges to the unique solution x =A“ for initial vector value x^ e R” if and 


only if the spectral radius p(/M™"N)<1 , where T - MN is called the iteration matrix. 


Here, we introduce some stationary iterative methods, shortly. 


Com. Alg. Num. Dim. 1(1) (2022) 1-24 


Edalatpanah | 


An experimental comparison of two preconditioned iterative methods to solve the elliptic partial differential equations 


2.1| Stationary Iterative Methods for SLE 
The Jacobi method ([5], [6]) 


For this method M=D and N=M-A. A single iteration step of the Jacobi method corresponds to the local 
solution of the linear system for a single variable. The resulting method is easy to implement, but often 


converges very slowly. 
The Gauss-Seidel method ([5], [6)] 


this is like the Jacobi method, except that it uses updated approximations as soon as they are available. In 
general, it will converge faster than the Jacobi method. For this method we have: M=D-L. 


The SOR method ([5]-[7]) 


This can be derived from the Gauss-Seidel method by introducing an extrapolation parameter w. For the 


optimal choice of w, the convergence speed can be increased substantially. The SOR method is defined by 
M 2 -(D - wL) And N=—((1-w)D + wu). 

So, we have Toog = (D -wL) [(1—w )D *wU ]. 

The SSOR method ([7]) 


Each iteration step of the Symmetric SOR (SSOR) method consists of two semi-iterations the first of 
which is a usual (forward) SOR iteration followed by a backward SOR iteration, namely an SOR where the 
roles of L and U have been interchanged. 


The AOR method ([8]) 


This method uses two parameters r (called relaxation parameter) and w (called extrapolation parameter). 


The AOR method is defined by M = -(D - rL) and N = - ((1 — w)D + (w — r)L) + wU). So, we have: 


T 
w 


=(D-rL)"[(1-w)D+(w-r)L+wU]. 


Tox 
The FIM method ([9]) 


Each iteration step of the FIM consists of 7 semi-iterations as follows: 


c P a+b) as 
X " =p [b+LX ^ +UX”], 
2 2 1 
(it+—) a (i+—) (i+—) 
n =p '[b+Lx ™+Ux 7] 
(4274) y, (+h) (+072) 
x ? =p [b+Lx ^ «Ux ^] 


n-1 . n-2 


: ‘ (i4——) 
x) =p [b+ Lx“? -Ul(1-w)x ^ +wx ^] .weR 


2.2| Non-Stationary Iterative Methods 


Non-stationary iterative methods differ from stationary methods in that the computations involve 
information that changes at each iteration step. This continuously updated information consists primarily 
of inner products of residuals or other vectors arising from the method. The best known non-stationary 
iterative methods for solving linear systems ate: 


The Conjugate Gradient method (CG) generates a sequence of conjugate vectors which are the residuals 
of iterates. These vectors are also the gradients of a quadratic functional, the minimization of which is 
equivalent to the solution of the linear systems. The CG method is extremely effective when the 
coefficient matrix is symmetric and positive definite [11], [12]. 


The MINRES method and the SYMMLQ method are alternatives to the CG method which are used if 
the coefficient matrix is symmetric but possibly indefinite. The SYMMLQ method generates the same 
solutions as the CG method if the coefficient matrix is symmetric and positive definite [5], [14]. 


The CGNE method and the CGNR method are specific CG methods for problems with non-symmetric 
and non-singular coefficient matrices. These methods are based on the fact that the A’A,AA’ and 
AA’ are always symmetric and positive definite. The CGNE method solves the system AA’Y =b for 
any y and then computes the solution x = A'Y . The CGNR method solves(A’A)x 2 b. for x where 


b=A'b . The convergence of these methods may be slow since the spectrum of A’A,AA’ and 
AA’ will be less favor able than the spectrum of A [5], [13]. 


The GMRES method computes a sequence of orthogonal vectors (as in the MINRES method), and 
combines them using a least-squares solve and update. However, unlike the MINRES method, it 
requires storing the whole sequence, so a large amount of storage is needed. This method is useful for 
general non-symmetric matrices [12], [14]. 


The BiCG method generates two sequences of vectors: one based on a system with the original matrix 
A and one A’, which are made mutually orthogonal, or bi-orthogonal. The BiCG method is useful when 
the matrix is non-symmetric and non-singular [12], [13], [15]. 


The QMR method applies a least-squares solve and update to the BICG residuals, thereby smoothing 
out the irregular convergence behavior of the BICG method. It also largely avoids the breakdowns that 
can occur in the BICG method [13], [15]. 


The CGS method is a variant of the BICG method that applies the updating operations of both 
sequences to the same vectors. An advantage is that this method does not need the multiplications by 


AT, but on the other hand convergence may be much more irregular than for the BICG method [13], 
[14], [16]. 


The BiCGSTAB method is a variant of the BICG method, like the CGS method. The difference is that 
the BICGSTAB method uses different updates for the sequence corresponding to AT in order to obtain 
smoother convergence than the CGS method [13], [14], [17]. 
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3 | Preconditioning Methods for Systems of Linear Equations 


A matrix P is called a preconditioner if the matrix P enables a transformation of the linear system into an 
equivalent system (with the same solution) which has more favorable spectral properties. For Eg. (1) 
preconditioning, transforms the system to 


Furthermore, it can be transformed to 


AFy=b, x=Fy, (7) 


where P and F are linear operators, called left and right preconditioners respectively. And for example, for 
stationary iterative methods we have: 


x") =MÍN x? - M?Pb, i=0,1,.... 
P P P 


where PA=Mp-Np and Mp is nonsingular. 


Also 


yo? = M;N,y” + Mjb, i- 0,1,.. 


where -AF-Mr-Nr and Mr is nonsingular. 
The majority of preconditioners fall in the first category; i.e. Eq. (6). 


The purpose of preconditioning is to change the matrix of the system, in order to accelerate the 
convergence of iterative solvers. 


We note that applying a preconditioner involves extra cost and most preconditioners involve an amount 
of computational work proportional to the number n of variables in their application. They thus multiply 
the amount of computational work per iteration by a constant factor. On the other hand, the number of 
iterations is usually only improved by a constant, i.e., the saving in computational works is independent of 
the matrix size. Therefore, there is a trade-off between the cost of constructing and applying the 
preconditioner and the gain of increased convergence speed. 


To improve the convergence rate of a basic iterative method, various models of preconditioning systems 
have been proposed. Here, we introduce some models of preconditioning, shortly. 


Jacobi Preconditioning ([18],[19]) 
The Jacobi or point-preconditioner consists of just the diagonal of the matrix A: 
P := diag( a11, a22, ..., Ann): 
Theoretically, there is no need for any extra storage, but in practice storage is allocated for the reciprocals 


of the diagonal elements in order to avoid repeated (un-necessary) division operations. 
If the index set J:={7,...,7}is partitioned into disjoint subsets J;, a block version can derived: 


ofai ^ if if Ik 
Pä = |o 


else: 
The preconditioner is now a block-diagonal matrix; this is the Jacobi block-preconditioning. 
Jacobi preconditioners need very little storage, and they are easy to implement, even on parallel 
computers. On the other hand, more sophisticated preconditioners usually yield much better 


improvements in the rate of convergence. 


SSOR preconditioning ([5], [20]) 
If the symmetric matrix A is decomposed as A = D + L LT, then the SSOR matrix is defined as: 
P:=(D+L)D -1(D + L5 


ot, parameterized by w, 


1 [1 MER T 
s, zl a a: 
2- w \w w w 
Using the optimal value of the relaxation factor w reduces the number of iterations, but in practice, this 
optimal value is prohibitively expensive to compute. As P is given in factored form, this method shares 
many of the properties of preconditioning methods based on incomplete factorizations. Since the 
factorization of the SSOR preconditioner is given a priori, there is no danger of a breakdown in the 
construction phase. 


Note. Above preconditioners called also preconditioners based on relaxation technique. For example, if 
the preconditioner P be as P=D, then this preconditioner is called Jacobi. Also if P = -(D - wL) then 


we have SOR preconditioner where for w = 7, we have Gauss-Seidel preconditioner and so on. 
Incomplete factorization [13], [20]- [23] 


Originally, preconditioners were based on direct solution methods in which part of the computation in 
skipped. This leads to the notion of incomplete LU (or ILU) factorization. 

The construction of an incomplete factorization may break down (due to a division by zero pivot 
element). The existence of an incomplete factorization is guaranteed for many factorization strategies, if 
the original matrix A has certain properties [23], [24]. 

An important consideration for incomplete factorization preconditioners is the cost of the factorization 
process. These costs bear fruit only if the iterative method without preconditioning converges very slowly, 
ot if the same matrix P can be used for solving several linear systems, as can be done, for instance, in 


Newton’s method for large nonlinear systems with a sparse Jacobian matrix. 


ADI preconditioning 

In 1994, Starke proposed the use of the Alternating Direction Implicit (ADI) method [25] as a 
preconditioning technique for the Krylov subspace methods for non-symmetric linear systems. Let A = 
A, + Az,where A4and A; are nonsingular. The alternative direction implicit method for solving the linear 
system Ax = bis in the following form: 


(A y+ riu. i =b- (A> = riDu;, 
2 
(A 2t roDu;4 =b- (A i= ralu, 1. 
2 


The ADI preconditioner is as P := (A, + 1% 1)(Az + r4I) . Parameters 11,172 are acceleration parameters. 
Varga [6] and Young [7] presented the optimum value for ri, r2. 


preconditioners of (I+S)-type 


In the literature, various authors have suggested different models of ([+S)-/pe preconditioner for the 
above-mentioned problem. In 1987 Milaszewicz [26] presented the preconditioner I+$"), where the 
elements of the first column below the diagonal of A eliminate. Gunawardena, Jain and Snyder 
considered [27] a modification of Jacobi and Gauss-Seidel methods and reported in 1991 that the 
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convergence rate of Gauss-Seidel method using the following preconditioning matrix is superior to that 
of the standard Gauss-Seidel method . Where, 


P=I+S. 


And, 


-a., for j=i+1,i=1,2,..,n-1 
sels Jj 0, for otherwise. 


Kohno et al. [28] proposed an extended modification of Jacobi and Gauss-Seidel methods. Their 


preconditioner is (J +s), where 


=Q a; j=i+1 Usa <1 


ES 0, otherwise. 


They also presented some numerical investigation for the choice of the optimal parameter. In [29] Usui et 
al. proposed to adopt 


P=I+U,(or I+L), 


as the preconditioned matrix, where U (L) is strictly upper (strictly lower) triangular of matrix A. They 
obtained excellent convergence rate compared with that by other iterative methods. Kotakemori et al. in 
[30] used 


P=I+S_, 
max 
where Smax is 


for 121,2,.,n-1,]»1 


Sa (s; ) - ^0. © for otherwise, 
and, 

V, «min je imax for i<n. 
Harano and Niki [31] considered the preconditioner 

P=1+(1+y)(L+U) 


where y is small positive number. 


Hadjidimos et al. [32] extended, generalized and compared the previous works. Wang and Song [33] 
presented a general form of the preconditioners for nonsingular M-matrices. Saberi Najafi and Edalatpanah 
in [34] established 


PzI-*S,.;; 


where for i=1:n-—1 and j»i: 


0, if jeQ, 
Pio Cà for otherwise, 
and, 
Q, -4j [jj«i&l, |= minja,, for i«n- 1. 


Furthermore, some other researchers have considered different models in the literature [35]- [44]. 


3 | Results 


In this section, we apply some iterative methods for solving elliptic PDE with Finite Differences Methods 


(FDM). Here, we establish our research about this topic. 


Consider the following elliptic PDE (advection-diffusion equation): 


ou ou 
(2-5) tent 


With the following bonditions; 
Dirichlet on x=0, y=0, either Neumann at yz Y and x=X or Dirichlet on y=Y and x=X, 


where 
e>0;xe[0,X] and y e| 0,Y | a(x,y)=1+x? ;b(x;y - Xe", 


and, 


x x y x 2 x y 2 
LXV) = € (rex +e - E NN J) pag x (Xe — we) E e&t QUT xev Hye). 


2 x2 


This equation has the following analytical solution: 
u(x,y) =o" (1 =e Jy 


Based on FDM, we solve the equation after the following differencing: 


I. Central differencing. 
II. Upwind differencing. 
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To see the convergence behavior of above FDMs by two iterative methods, we select Saberi Najafi and 
Edalatpanah’s method (FIM [5]) as stationary iterative method and a good preconditioned nons-tationary 
iterative method called the preconditioned (Incomplete LU Decomposition) BICGSTAB method. And 
finally, we show: 


— Plot the residual reduction. 
— Plot the error map. 
— Estimate the truncation error. 


— Plot the smoothing factor of the iteration. 


For the above problem, two types of discretization schemes (centered and hybrid) have been compared 
for the comparison of the solution. 


For both schemes, we choice: 


= iis Uiaj A TH E O(A x), 

ay echoes aga) 
And for hybrid Scheme: 

^q = e O(a, 

Adj = + O(y) 


Then by above equation, we get: 


[e Ax? - byAx* Ay] Ui 4 + [E Ay? =a, AxAy^ | Uia; [ 2€ Ay? -2e Ax? + rag Ax Av" 
t bj Ax?Ay | Ui; +[e Ay? ] Uia +[[e Az ] Uija = fi Ax? Ay’, 


where, 
aj =1 + x?, by =Xe % 


Furthermote, for centered Scheme: 


9u,.. | Uiaj-Uiij 
E. jij = Lum m O(A x) 
ðu. | Uiji-Uiji 
ax id 7— 34, — + O(Ay) 


So we obtain, 


[e Ax? = 0.50b; Ax? Ay ] Uiia + [ € Ay? = 0.50a;j AxAy? ] Ui; + [ -2€ Ay? - 2€ Ax? ] 
Ui + [ E€ Ay? + 0.50a;j AxAy? ] Ui; + [[ E€ Ax?+ 0.50b; A x? Ay] Ui; i= fs Ax? Ay?. 


Now, we present the results: 


If we choose € 2 4, X =Y = 4, Tolerance=1e-20 and number of mesh size on both direction=20. Then 
with FIM 5-step with k=-3.7, we have: 


For dirichlet boundary condition 


Figs. 1-2, show the numerical solution of the mentioned equation with two types of discretization 
schemes. 


Numerical solution using FIM-5-step 


Numerical solution using FIM-5-step 


Fig. 1. Solution with upwind (hybrid) scheme. Fig. 2. Solution with centered scheme. 


Figs. 2-4, show the error plot for the solution of this problem by FIM and two types of discretization 
schemes. 


Error Plot for the solution using FIM-5-step 


Error Plot for the solution using FIM-5-step 


Error(x, y) 


Error(x,y) 


Fig. 3. Error plot with upwind (hybrid) scheme. 


Fig. 4. Error plot with centered scheme. 


For neumann boundary condition 


Figs. 5-6, show the numerical solution of the mentioned equation with two types of discretization 
schemes. 
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Numerical solution using FIM-5-step 


Numerical solution using FIM-5-step 


Fig. 5. Solution with upwind (hybrid) schem. Fig. 6. Solution with centered scheme. 


Figs. 7-8, show the error plot for the solution of this problem by FIM and two types of discretization 
schemes. 


Error Plot for the solution using FIM-5-step 


Error Plot for the solution using FIM-5-step 


Fig. 7. Solution with centered scheme. 


Fig. 8. Solution with centered scheme. 


In Table 1 , the convergence behaviors of the FIM and a good non-stationary preconditioned method called 
ILU preconditioned BiCGSTAB method (Pre- BiCGSTAB) are compared by CPU time, number of 
iterations (Iter) and the maximum error. 


Table 1. Comparison of the results between two different iterative methods. 


Dirichlet Boundary Condition Neumann Boundary Condition 
Hybrid Scheme Centered Hybrid Scheme Centered Scheme 
Scheme 
FIM 23 24 128 136 
Pre- 35 39 47 55 
BiCGSTAB 
Maximum FIM 0.0073 1.1829e-04 0.0113 0.1252 
E 
ind Pre- 0.0073 1.1404e-04 0.0110 0.1256 
BiCGSTAB 
CPU time FIM 0.1643 0.1688 0.0118 0.3668 
Pre- 2.9106 3.2193 3.8977 4.5627 


BiCGSTAB 


From the result we can see that: 


As mentioned in the above table, iteration required for the Neumann case is more as compared to its 
Dirichlet counterpart. 

For the Dirichlet boundary condition, the number of iteration and CPU time for the hybrid scheme is 
less than Centered scheme. However, the maximum error for the centered scheme is less. 

For the Neumann condition, the number of iterations, maximum error and CPU time for the hybrid 
scheme is less than centered scheme. 

FIM method is superior to the ILU preconditioned BiCGSTAB method. 

For low value of epsilon, both methods don't work. The matrix may be become ill-conditioned. To 
show the convergence behaviots of these methods, here we test the problem with different parameters 
of e for Dirichlet boundary condition. 


Table 2. Results by both methods for epsilon=0.0001. 


e—0.0001 Hybrid Scheme Centered Scheme 


FIM 


Pre- 
BiCGSTAB 


Solution oscillates and 
does not converge after 
1000 iteration 


Error Plot for the solution using FIM-S-step 


Iter =2, Cpu time = 0.6075 
Maximum_Error = 0.0406 


Error Plot for the solution using ILU(D)preconditioned BICGSTAB 

l Solution oscillates and 
does not converge after 
1000 iteration 


Iter =30, Cpu time =2.6066 
Maximum_ Error =0.0406 
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Table 3. Results by the FIM method for epsilon=0.001. 


e=0.01 Hybrid Scheme Centered Scheme 


Error Plot for the solution using FIM-5-step. 


Solution oscillates and does not 
converge after 1000 iteration 


FIM 
Iter =5, Cpu time = 0.13817 
Maximum, Error = 0.04148 
Table 4. Results by the Pre-BiCGSTAB method for epsilon=0.001. 
e=0.01 Hybrid Scheme Centered Scheme 
Error Plot for the solution using ILU(0)preconditioned BICGSTAB 
i Solution oscillates and does not 
converge after 1000 iteration 
Pre- 
BiCGSTAB 
Iter =29, Cpu time = 2.47187 
Maximum_ Error = 0.041448 
Table 5. Results by both methods for epsilon=0.1. 
e=0.1 Hybrid Scheme Centered Scheme 
FIM Solution oscillates and does not Solution oscillates and does not 


converge converge 


Error Plot for the solution using ILU(D)preconditioned BiCGSTAB 


3 


Solution oscillates and does not 


Pre- converge 
BiCGSTAB 

Iter =91, Cpu time = 7.636989 

Maximum, Error — 0.0018694 

Table 6. Results by FIM method for epsilon=1. 

e=1 Hybrid Scheme Centered Scheme 

Error Plot for the solution using FIM-5-step 

Solution oscillates and does not MEE 
Error Plot x10 
converge 
[ES 
54 
ü 
Iter 219, Cpu time = 1.653104 
MTS Maximum, Error — 4.57593229e-04 
Residual , Residual Reduction Plot 
Reduction a a ME 07 
Plot “| 
ash 
07t 


Norm of Residue 
[= 
an 
T 


10 12 14 15 18 20 
No of Iteration 
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Table 7. Error plot of the Pre-BiCGSTAB method for epsilon=1. 


e-1 Hybrid Scheme Centered Scheme 


Error Plot for the solution using ILU(D)preconditioned BICGSTAB 


Error Plot for the solution using ILU(D)preconditioned BICGSTAB 


Error E ot 
Plot 1 
Iter 2228, Cpu time= 18.9332 Iter 219, Cpu time= 1.6502 
Maximum, Error— 0.0235 Maximum, Error- 4.5759e-04 
Table 8. Residual plot of the Pre-BiCGSTAB method for epsilon=1. 
e=1 Hybrid Scheme Centered Scheme 
Residual Reduction Plot Residual Reduction Plot 
9000 r : ; 1 1 = : 
8000 - | ost J 
7000 4 I | 
Residual 6000 - 4 " d | 
Reduction E 3 8f ] 
$ 5000 | E 
Plot x s O57 ] 
g 4000} 1 5 oat J 
Z a0} d 03 J 
2000} J 02 J 
1000} E 21 | 
Se 09 3 4 6 68 10 1214 1 18 2 
0 50 100 150 200 250 No of Iteration 
No of Iteration 
Table 9. Results by the FIM method for epsilon=2. 
£= Hybrid Scheme Centered Scheme 
Error Plot for the solution using FIM-5-step 
Error Plot for the solution using FIM-5-step 4 
Error Plot m 
Residual 
Reduction 
Plot Iter =93, Cpu time = 0.301401397 Iter =18, Cpu time = 0.1635525769 


Maximum_Error = 0.0129400359 


Maximum, Error = 2.398448563e-04 


Residual Reduction Plot Residual Reduction Plot 
T T T T T T 


1 T T T 


09H J 
08H E 
07H E 
2 osp 4 $ 
E] d 
E 5 
= 05 J z 
5 gal 4 = 
z 
03- E 
02- 4 
01 E 
ü — o 2 4 6 8 10 12 14 16 18 
o i 2 3 40 5 6 70 8 90 100 No of Iteration 
No of Iteration 
Table 10. Results by the Pre-BiCGSTAB method for epsilon=2. 
e=2 Hybrid Scheme Centered Scheme 
Error Plot for the solution using ILU(D)preconditioned BiCGSTAB Error Plot for the solution using ILU(D)preconditioned BICGSTAB 
4 
x10 
Error Plot 
0.015 5 
_ 001 " 
oo Š 
4 
y 00 x 
Residual Iter =32, Cpu time = 2.72658957 Iter =28, Cpu time = 2.35891506 
Reduction Maximum, Error = 0.0128128 Maximum, Error = 2.82797714639e-04 
Plot Residual Reduction Plot Residual Reduction Plot 
35 : - - - - : , 
og E 
3r E 
08 4 
25} E 07 J 
2 2 06 J 
$ 2 | i 
= = os J 
E 15F 4 E oat 4 
5 2 
= 
ost 4 
1b 4 
0.2} 4 
05r 4 01 J 
0 A 
o " i um n] 5 10 15 20 25 30 
0 5 10 15 20 25 30 35 No of Iteration 


No of Iteration 


By above results we can conclude that: 


II. 


In both methods, for the large value of e, the convergence of the solution is much more ensured. 
The etror is less for large e. The stiffness matrix is well conditioned for large values of e. The stiffness 
matrix is much mote diagonally dominant for the case of higher values of epsilon. 

In both methods, for very low values of e, the centered scheme does not converge. In general, the 
hybrid scheme converged with very less iteration and the centered scheme converged but with large 
number of iteration and with high oscillation before convergence. 
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HI. Convergence area of the ILU preconditioned BiCGSTAB method is superior to the FIM method. 
Howevet, the speed of convergence of the FIM is the best. 
IV. In both methods, for € € (.01,0.9), the hybrid scheme does not converge. However, the centered 
scheme converged. 
V. For higher value of e, both the scheme converges. However, the error in case of centered scheme is 
much less as compared to the hybrid scheme. 


Next, bye =4,X =Y = 4, Tolerance-1e-20 and number of mesh size on both direction=20, we show 


the residual reduction plots for some iterative methods and hybrid discritization scheme (see Tables 11-14). 


Table 11. Residual plot of the ILU method [12]. 
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Table 12. Residual plot of the BICGSTAB method. 
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Table 13. Residual plot of the Pre-BiCGSTAB method. 
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Table 14. Residual plot of the FIM method with k=-5.5. 
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Next, we compare the methods based on amplification factor plot (flow parameter is same as before). 


kl 
d 


k 
ul 


The amplification factor is expected to converge to a value less than 1(see Tables 15-18). 


The amplification factor is computed as: 


, Where the index & denotes the number of iterations. 
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Table 12. Amplification factor plot of the ILU method [12]. 
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Table 18. Amplification factor plot of the FIM method with k=-5.5. 
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Table 19, shows a total comparision between two models. 


Table 14. Brief statistics. 


ILU[7] BiCGSTAB Pre-BiCGSATB FIM_5-step(K=-3.7) 
No. of Dirichlet NC 91 35 23 
Iteration 
Required | Neumann | NC 118 47 128 
(with cpu time= (with cpu (with cpu time=0.0118) 
7.5998) time=3.8977) 
Maximum | Dirichlet NC 0.0073 0.0073 0.0073 
Error 
Neumann NC 0.0110 0.0110 0.0113 
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Fig. 1.CPU time for non-preconditioned BiCGSTAB 
(Hybrid scheme). 


Fig. 2 .CPU time for preconditioned BICGSTAB 


(Hybrid scheme). 
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In all of above results the numbet of mesh size on both direction is 20. Here, we study the results based 
on different number of mesh size. In Fzgs. 9-12 we show CPU time of different methods with different cell 
number for hybrid discritization scheme and Dirichlet boundary condition. In Figs. 13-16 we show CPU 
time of different methods with different cell number for centered discritization scheme and Dirichlet 
boundary condition. 
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Fig. 3. CPU time for FIM_5-step (Hybrid scheme). Fig. 4. CPU time for FIM_7 step(k=-5.5) (Hybrid 
scheme). 
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Fig. 5.CPU time for non-preconditioned 


Fig. 6 .CPU time for preconditioned BiCGSTAB 
BiCGSTAB (Centered scheme). 


(Centered scheme). 
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. 7. CPU time for FIM. 5-step (k=-3.7) (Centered Fig. 8. CPU time for FIM 7 step(k=-5.5) (Centered 
scheme). scheme). 


From that above figures we an see that: 


I. The time required for case of FIM. 7-step is less as compared to all other cases. 
IL The maximum time required is for the case of BICGSTAB. 
III. Preconditioned BiCGSTAB takes less time as compared to the No preconditioned BICGSTAB. 
IV. FIM. 7-step takes less time as compared to the FIM_5step. 
V. Generally, hybrid discritization scheme takes less time as compared to the centered discritization scheme. 


4 | Conclusion 


In this paper, we investigated the comparsion of two iterative methods for solving the elliptic partial 
differential equations. Our expriments can be explained as follow: 


I. Among all the methods discussed, FIM method plays well as per as CPU time and convergence and 
other features ate considered. 


IL | Large number of cells gives better result. But it takes more time to converge. 

IIl. Preconditioned used also plays nice role in the CPU time consumption and the converges of the 
solution is considered. 

IV. The value of e plays a crucial role as depending on the e, the contribution of the diffusion term will 
determine. 
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